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ABSTRACT 



Infrared and submillimetre observations of nearby Vega-like stars have revealed a 
number of clumpy, asymmetric dust debris disks. Previous studies using semi-analytical 
and numerical methods have suggested planetary companions of various mass as the 
likely cause of most examples of disk asymmetry. In this paper, we modify an N-body 
fT) \ symplectic gravitational integrator to include radiation terms and conduct medium- 

resolution parameter searches to identify likely planetary candidates in observed Vega- 
© 1 like systems. We also present high resolution models of Vega and e Eridani, comparing 

^ ■ our results to those of previous authors, and a new model for Fomalhaut. 

O , Subject headings: circumstellar matter - methods: N-body simulations - planetary 

q . systems - stars: individual (Vega, e Eridani, Fomalhaut) 



1. Introduction 

Dust disks composed of particles in the micron to centimetre range are known to exist around 
many stars that exhibit an infrared excess, such as Vega, (3 Pictoris, e Eridani and others (see e.g. 
Zuckerman 2001). Particles of this size are affected significantly by solar radiation and corpuscular 
(solar wind) forces, with Poynting-Robertson (PR) and solar wind drag resulting from the moving 
particle's absorption and reemission of solar energy. Over time, these dissipative forces caused by 
radiation and solar wind act to reduce a dust particle's semi-major axis and eccentricity, and the 
particle eventually spirals into the central star (Burns et al. 1979). Since the age of most of these 
systems exceeds the expected lifetime of a small dust particle, the particles must be continually 
replenished. The suspected mechanism is a collisional cascade involving primordial planetesimal 
and smaller sized bodies (Backman & Paresce 1993; Wyatt & Dent 2002) as is believed to occur 
in the Edgeworth-Kuiper Belt (e.g. Backman et al. 1995; Landgraf et al. 2002). The destructive 
process which leads to dust formation has led these disks to be known as 'debris disks'. 



Of the debris disks near enough to Earth to be spatially resolved, most show a significant degree 
of asymmetry (e.g. Vega: Holland et al. 1998; e Eridani: Greaves et al. 1998). The commonly 



- 2 - 



held view is that gravitational interactions with an unseen embedded planet are responsible for 
most cases of disk structure (Ozernoy et al. 2000; Quillen & Thorndike 2002; Wilner et al. 2002). 
Alternative possible explanations have been considered for specific systems, such as disturbance by 
a passing star (e.g. HD141569: Clampin et al. 2003), a recent large cometary/planetesimal collision 
(e.g. Fomalhaut: Wyatt & Dent 2002) and contamination by a background feature such as a distant 
galaxy (e.g. Vega: Wilner et al. 2002). 

A planetary companion has several effects on an initially uniform debris disk. Particles near 
the planet are often ejected due to a close encounter - an effect which is enhanced as the planet's 
mass increases. Since particles interior to the planet spiral increasingly rapidly into the central star, 
an interior cleared zone is usually present (Liou & Zook 1999). However, some particles spiralling 
in towards the star may become trapped in a series of Mean Motion Resonances (MMRs) with the 
planet. An MMR exists when the ratio of orbital periods for two objects orbiting a central mass is 
equal to m : n, where m and n are integers. The order of an MMR is given by the quantity \m — n\. 
Earth has a dust ring formed by dust particles in MMRs (Dermott et al. 1994), and Neptune is 
suspected to do the same to Kuiper Belt dust (Liou & Zook 1999). 

Whilst locked into an MMR with a planet, the angular momentum of the particle lost to 
Poynting- Robertson and solar wind drag is balanced by the resonant forcing of the planet's gravi- 
tation. A particle trapped in an MMR can extend the particle's lifetime to many times the value 
predicted by Poynting-Robertson and solar wind drag alone. The resultant 'pile-up' of particles at 
discrete semi-major axes tends to cause prominent radial and angular structure in the debris disk 
(Kuchner & Holman 2003). 

In multiple planet systems, it is the outermost planet that generally dominates the structure 
of the resulting debris disk. Neptune is suspected to play this role in interactions with Kuiper Belt 
dust (Liou & Zook 1999). The interior planets, however, can have significant effects on particles 
which 'leak past' the outermost planet. In the solar system, for example, Jupiter and Saturn 
eject many dust particles which drift interior to Neptune (Liou & Zook 1999). For simplicity, we 
investigate only single planet systems and assume that the effects of any interior planets which may 
be present is negligible outside the orbit of the outermost planet. 

In this paper, we numerically simulate a large number of debris disk systems, creating a 
synthetic debris disk catalogue. Our numerical method is described in Section 2, and the testing 
of the code is detailed in Section 3. We present a small sample of the results from the synthetic 
catalogue in Section 4, and use the catalogue to select planetary configurations which could be 
responsible for the observed systems of Vega, e Eridani and Fomalhaut. The selected systems are 
then simulated in higher detail in Section 5, and the results compared to observational data. Our 
conclusions are presented in Section 6. 
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2. Numerical Method 

We model the evolution of debris disks using an N-body integrator that includes the effects 
of radiation pressure, Poynting-Robertson and solar wind drag. In this section we discuss the 
modifications made to a symplectic integrator to account for the additional forces, as well as our 
numerical techniques for creating debris disks and recording particle positions. 

2.1. The RMVS3 Integrator 

Since the general problem of a particle moving under the influence of gravitational and radia- 
tion forces is nonlinear, it is necessary to numerically integrate the orbits of all particles. Commonly 
used integrators include Runge-Kutter, Burlisch-Stoer and Mixed Variable Symplectic (MVS) al- 
gorithms. The MVS integrator developed by Wisdom & Holman (1991) is normally the fastest, 
but has the disadvantage of being unable to follow close encounters between particles and planets, 
which are crucial to debris disk evolution. This drawback has been overcome with the use of in- 
tegrators such as RMVS3 (Levison & Duncan 1994) and SyMBA (Duncan et al. 1998), which are 
based on the Wisdom & Holman (1991) scheme but can integrate close encounters. In this work we 
have modified RMVS3 to include radiation and solar wind forces which affect debris disk particles. 
Such an approach was used by Moro-Martin & Malhotra (2002), who modify a variant of SyMBA. 
(Note that SyMBA and RMVS3 handle close encounters between particles and planets in different 
ways.) 

RMVS3 assumes that test particles are both massless and collisionless - the problem of dust 
collisions is discussed in Section 4.1. The code requires units such that the gravitational constant 
G = 1 and so we have chosen distance, time and mass units to be 1 AU, 1 year and M /(47r 2 ) 
respectively. Throughout this paper we shall use a for the semi-major axis, e for eccentricity, and i 
for inclination, with the subscripts £p, pb and pi referring to a test particle, parent body and planet 
respectively. 

As developed in Wisdom & Holman (1991) and described in Levison & Duncan (1994), RMVS3 
expands the Hamiltonian of a test particle into two integrable components given by 

H = H kep + Hdist ? (1) 

where H kep represents the Keplerian motion around the central star, and Hdist represents the 
perturbances on a test particle resulting from the planet (s). Using a timestep At, the second order 
approximation implemented in RMVS3 consists of applying the disturbance Hamiltonian for At/2, 
applying the Keplerian Hamiltonian for At, and then applying the disturbance Hamiltonian again 
for At/2. The approximation is accurate as long as the planetary disturbances are small compared 
to the Keplerian evolution. This assumption of small disturbances does not hold when a particle 
suffers a close encounter with a planet. In this situation, since the gravitational effect of the planet 
on the particle exceeds that of the star, the roles of planet and star are reversed and the Keplerian 



-4- 



portion of the Hamiltonian is taken to occur around the planet and the star is relegated to the 
role of disturbance Levison & Duncan (1994). Although a smaller timestep is used during a close 
encounter, the duration of an encounter is typically small enough that the simulation proceeds 
rapidly. In the intermediate region where planetary and stellar forces are comparable, the timestep 
is also reduced but the integration remains heliocentric. 



2.2. Addition of Radiation and Solar Wind 

The magnitude of the radiation force on a particle is given by 

lFradl = ' {2) 

where L is stellar luminosity, A is the particle cross-sectional area, Qpr is the radiation pressure 
coefficient, c is the speed of light and r is the heliocentric distance. It is standard practice to 
describe the strength of the radiation force on a given particle as a fraction of the gravitational 
force on a particle, such that 

Frad — /3Fg rav . (3) 

The constant (3 is depends on the particle size and composition, as well as the stellar mass and 
luminosity, but is independent of the particle's heliocentric distance. The magnitude of Poynting- 
Robertson drag is proportional to (3. The ratio of solar wind drag to Poynting-Robertson drag is 
dependent on particle size, but is relatively constant over the range of particles usually considered 
(s > 0.5/xm, Burns et al. 1979), so the combined drag effects of PR and solar wind can be expressed 
using a single parameter given by 

Psw =f3(l + sw) , (4) 

where sw is the ratio of solar wind drag to Poynting-Robertson drag (Burns et al. 1979). Following 
Moro-Martin & Malhotra (2002), we use a constant value of sw = 0.35 in all our simulations. 

The acceleration on a particle due to solar wind drag and radiation, ignoring terms of order 
(v/c) 2 and higher, is given by 

d 2 v 



dt* ~ Ff3rav 



n~ fisw / ~ . \ 

pr (v r r + v) 

c 



(5) 



The first term in this expression is the result of the outwardly-directed radiation pressure, which 
causes the net acceleration in the radial direction to be reduced to (1 — f3)F grav . The remaining two 
terms (which are smaller by a factor of v/c) represent the Poynting-Robertson and solar wind drag. 
Poynting-Robertson drag results from the particle's motion relative to the radiation source, which 
causes a component of the radiation force to oppose the particle's motion. During normal particle 
movement, the first term is accounted for in the Keplerian Hamiltonian by reducing the central 
object mass to the 'apparent' value, while the other terms are added to the disturbance Hamilto- 
nian. During a close encounter with a planet, the entire expression is added to the disturbance 
Hamiltonian. 
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2.3. Particle Initialisation 



All simulations commence with a central star, an orbiting planet and a number of test particles 
distributed around the star with a fixed value of /?. The star is fixed at the origin and has mass 
represented by M*. The orbiting planet is described by M p i, a p i and e p i, representing the planet's 
mass, semi-major axis and eccentricity respectively. 

Test particles are assumed to be released from larger 'parent bodies' (which are unaffected by 
radiation pressure) at the beginning of the simulation. This is analogous to the creation of dust 
through the collision of larger rocky bodies in a debris disk. The concept of a parent body is used 
only to set initial orbital parameters for test particles; they play no part in the simulation after 
initialisation. For each simulation, the initial distribution of parent bodies is specified by a range 
of semi-major axis, eccentricity and inclination values, given by a^, e p & and respectively. The 
arguments of perihelion, longitudes of ascending nodes and mean anomalies for the parent bodies 
are allocated randomly. A test particle is then created from each parent body, with semi-major axis, 
eccentricity and inclination referred to by a tv ,e tv and i tp respectively. Since the test particles are 
affected by radiation pressure while the parent bodies are not, the test particles will have different 
orbital parameters from their parent bodies. In particular, the test particles will have increased 
semi-major axes and changed eccentricities, given by 



1-/3 



etp 



= 1 



2a pb (5/r' 

(l-2a pb p/r)(l-e 2 pb ) 
(1 - (5f 



1/2 



(6) 
(7) 



2.4. Particle Recording 

Simulating a sufficient number of test particles to resolve the disk structure at all times is com- 
putationally expensive and an alternative method, used by Liou & Zook (1999) and Moro-Martin & 
Malhotra (2002), is to simulate a smaller number of particles and regularly record their positions, 
which can then be summed to produce distribution maps. A simulation is terminated after a fixed 
time or when all particles have been destroyed through either ejection via close encounters with 
planets or accretion onto the central star. In this way, long lived particles (trapped in resonances) 
contribute more to the final dust distribution. If we assume the dust distribution is ergodic and that 
the total dust mass is constant (i.e. the dust production rate equals the loss rate), this approach 
should give an accurate representation of the actual 'steady-state' disk structure. 

In previous studies (e.g. Moro-Martin & Malhotra 2002) the particle locations have been trans- 
formed into a reference frame which rotates with the outermost planet, but this approach is not 
appropriate for planets with eccentric orbits. For this reason, we instead ensure that the particle 



-6- 



recording takes place with the planet in specific orbital phases. In this way a number of 'snap- 
shots' of the system are taken with the planet in varying locations, allowing the effect of planetary 
phase to be studied (Wilner et al. 2002; Kuchner & Holman 2003). For example, many simulations 
utilised a planet with orbital period 300 years and a particle recording time of 1000 years, resulting 
in particle distributions for three different planetary phases. 

It should be noted that there are some limitations to this procedure of simulating a debris disk 
using a small number of test particles. As discussed by Moro-Martin & Malhotra (2002), because 
both the probability of being trapped in a particular resonance and the time spent in a resonance 
depend sensitively on initial conditions, using a limited number of particles may over-exaggerate the 
importance of a particular resonance. We have tried to overcome this problem by using an order of 
magnitude more particles than previous simulations (e.g. Quillen & Thorndike 2002; Moro-Martin 
& Malhotra 2002) in high resolution simulations (see following sections). An additional problem 
with this procedure is that all of the test particles are released from their parent body at the same 
time, with a specific planetary phase. In reality, dust particles would be released continuously, at 
random planetary phases. This effect, however, is not likely to have a significant impact, since 
the test particles are typically released a long way from the planet and are distributed randomly 
around the star. 

Setting a fixed maximum simulation time may also be problematic. For the small particle 
number approach to be accurate, the integration should be continued until all test particles have 
been destroyed. The termination of a simulation at a predetermined time means that some particles 
may still be active at the end of the simulation. One could argue, however, that extremely long- 
lived particles are unlikely to exist due to destruction via particle collisions. Since RMVS3 is a 
collisionless code, the destruction of particles in this way is not accounted for. We investigate these 
effects in more detail in Section 4.1. 

3. Test Calculations 

There is no general solution for the motion of a particle under the influence of radiation and 
gravitational forces in a system comprising a star and one or more planets. However, analytic 
solutions and analytic approximations exist for particle motion in the two-body (Wyatt & Whipple 
1950) and the circularly restricted three-body (Liou & Zook 1997) problems respectively. These 
specific cases can be used to check the accuracy of the modified RMVS3 code. 

To ensure that the addition of radiation and solar wind forces did not effect the normal running 
of the code, we ran a series of tests in which (3 was set to zero, and compared the results to identical 
tests run using the original RMVS3. The results were indistinguishable in all cases. 
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3.1. The 2-body problem 

The time rate of change of a test particle's orbital elements in the radiation 
problem are given in Moro-Martin & Malhotra (2002) (following the work of 
(1950); Burns et al. (1979)), and are equal to 

fda\ _ _(l + sw)pM* (2 + 3e 2 ) 

\dt) ~ c a (l- e 2) 3 / 2 ' 

fde\ _ h{l + sw)(3M+ e 

\Jt) ~ 2~c a 2(l-e 2 ) 1/2 ' 

(!) - 

where is the mass of the central star, and c is the speed of light. 

To test the code against the two-body solutions, we ran a simulation with a single /? = 0.15 test 
particle placed in orbit around a solar-mass star, released from a parent body with a p b = 250 AU, 
e pb — 0.8 and i p b = 7.6°. The particle's argument of perihelion, longitude of ascending node and 
mean anomaly were selected randomly. The analytical and numerical results plotted in Figure 1 
are in excellent agreement. 



-modified two body 
Wyatt & Whipple 

(8) 
(9) 
(10) 



3.2. The restricted 3-body problem 

The addition of a planet greatly complicates the motion of the test particle, but Liou & Zook 
(1997) have derived expressions for the time evolution of a particle's orbital elements whilst in an 
MMR with a zero eccentricity planet. These expressions are a second order approximation, valid 
only for low test particle eccentricities and inclinations, and low order resonances, and are given by 

e 2 = (el-(K-l)/3))exp(3At/K) + (K-l)/3, (11) 
i = i exp(-At/4), (12) 

where A = (2(3 SW M+) /(a 2 c), K — m/n, eo and zq represent the particle's initial eccentricity and 
inclination, and m and n are the integers specifying the resonance. It should be noted that the 
location of resonances is shifted from their normal location due to the presence of the radiation 
force, which causes particles at a given semi-major axis to orbit more slowly due to the reduced 
effective gravitational force. The semi-major axis of a particle locked in an m : n resonance with a 
planet is given by 

a tp = a pl (1-/3) l ' z (m/nf 3 . (13) 

We start at time t = with a single particle placed exterior to the resonance in orbit around 
a one solar mass star and a planet at 30 AU. We follow the orbital evolution of the test particle 
as it passes through the MMR until it is ejected from the system, and compare the results with 
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equations 11 and 12. We show the results for two different j3 values. In Figure 2a, a single j3 — 0.05 
particle is trapped in a 3 : 2 resonance with a Neptune mass planet, while in Figure 2b a particle 
with /? = 0.2 is trapped in a 2 : 1 resonance with a Jupiter mass planet. The particles' orbital 
elements are well described until their eccentricity and inclination become too large and the analytic 
solution becomes invalid. The results shown in Figure 2 are only a small sample of the large number 
of tests carried out to ensure the accuracy of the code. 



4. A Synthetic Debris Disk Catalogue 

Previous studies by Kuchner & Holman (2003) have predicted the general effects of varying 
planetary mass and eccentricity on the dust distribution in a debris disk. We have generated a 
synthetic catalogue of debris disks which provides numerical agreement to the theoretical predictions 
and can be used as a guide when interpreting observed systems. Here, we present a sample of results 
obtained from modelling over 300 disks, showing the effects of planetary mass and eccentricity, 
dust composition (/?) and initial dust distribution on the disk structure. By generating simulated 
observations of these simulations, we show that even observations at a low resolution can distinguish 
between different planet and dust combinations. 

For the creation of our synthetic catalogue, we allowed five parameters to vary: M p i,e p i, /?, 
a p 5 and e p t,. All other parameters were fixed as follows: = 1M , 0° < i p b < 8°, n tp = 500, and 
planetary period = 300 years (yielding a p \ ~ 44.8 AU). Table 1 summarises the parameter space we 
have investigated. Of the 600 possible combinations of our five free parameters, approximately 300 
disks were simulated. Subsets of the parameter space which were fully explored include /? = 0.1, 
M p i = 0MM Jup and M pl = lM Jup . 



Fig. 1. — Time evolution of a single test particle orbiting a solar mass star subject to Poynting- 
Robertson and solar wind drag. Plotted are the evolution of (a) a tpi (b) e tp , and (c) i tp for a particle 
with (3 = 0.15 and sw = 0.35. The analytic solution (dashed line) and actual particle path (solid 
line) are coincident. 
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Fig. 2. — Evolution of a single test particle whilst in an MMR with a planet. The particle evolution 
is shown as a solid line, and the analytic solution is shown as a dashed line. The planetary semi- 
major axis is shown as a dotted line in the bottom two panels, (a) Particle with f3 = 0.05, in a 3:2 
resonance with a 0.05Mj up mass planet. The particle is trapped for about 10 7 years, (b) Particle 
with (3 = 0.2, in a 2:1 resonance with a Jupiter mass planet. The particle is trapped for about 5 
million years. In both cases sw = 0.35. 
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CLpb/ dpi Cpb 



M pl /M Jup 




P 






0.01 


0.0 


0.01 


[1.0, 1.4] (close) 


[0.0, 0.2] (low) 


0.05 


0.1 


0.1 


[1.4, 2.0] (mid) 


[0.1, 0.5] (medium) 


0.2 


0.3 


0.2 


[2.0, 3.5] (far) 




1.0 


0.5 


0.4 






3.0 


0.7 









Table 1: The range of parameter values used in generating our synthetic catalogue. Note that while 
this parameter space was not completely explored, approximately 300 of the 600 possible disks were 
simulated. Fully explored parameters (where the chosen parameter was simulated with all possible 
combinations of the other parameters) include j3 — 0.1, M p i = 0.05Mj up and M p i = lMj up . 
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Each simulation was run for 6 x 10 7 years, corresponding to 2 x 10 5 planetary orbits. Particle 
positions were recorded every 1000 years, allowing observation of the system in 3 different planetary 
phases. The maximum and minimum allowable heliocentric radii of the test particles were 700 and 
2 AU respectively. 

At the end of each simulation, the recorded particle positions (in the three orbital phases) 
were binned in the Xy-plane. The particle distribution in the Xy-plane was then plotted for each 
particular planetary phase. We also binned particle semi-major axes and plotted the semi-major 
axis occupancy versus semi-major axis, which shows the location and strength of MMRs with the 
planet. 

In order to produce synthetic observations of the systems, the particle distribution plots were 
weighted by estimated total disk mass and the dust emissivity at 850 /im, assuming the dust to 
be a perfect blackbody. This emissivity plot was then convolved with a two dimensional Gaussian 
with a full width half maximum (FWHM) of 40 AU, to simulate observations of the system with 
limited resolution. The synthetic observations can also be plotted from various disk viewing angles. 
This procedure assumes that the disk is optically thin, which is a reasonable assumption for disks 
containing little to no gas, as is the case for most Vega-like systems (Liseau 1999). 

The complete synthetic catalogue is available for viewing online 1 . Here we present a sample of 
results from our catalogue for Jupiter and Neptune mass planets with e p \ = 0.1 and 0.5 in Figures 3 
and 4. These simulations use /3 = 0.1, and test particles released from parent bodies with orbital 
elements < e p & < 0.2, 1.4 < a p b/a p i < 2.0 and < i p b < 8°. 

As can be seen in Figures 3 and 4, the mass and eccentricity of a planet dramatically alters 
the resulting dust distribution. The Jupiter-mass planet clears a central cavity, which is not as 
pronounced in the Neptune-mass case. In the low planetary eccentricity case, a relatively smooth 
dust ring is present external to the planet. When the planetary eccentricity is increased, the smooth 
ring vanishes, replaced by prominent arcs or clumps of dust emission. Similar trends were observed 
in the synthetic catalogue with planets of different masses. These results closely resemble the 
predictions made by Kuchner & Holman (2003) for the four cases of: (a) low mass, low eccentricity 
planets, (b) low mass, moderate eccentricity planets, (c) high mass, low eccentricity planets, and 
(d) high mass, moderate eccentricity planets. 

Figure 5 shows the resonance occupations for the e p \ = 0.1 cases. The effect of planetary mass 
is again highlighted, with a high mass planet trapping particles in a few strong resonances, mostly 
at higher semi-major axes, while a low mass planet has many weaker resonances, closer to the 
planet's semi-major axis. More accurate simulated observations for specific systems and telescopes 
are undertaken in Section 5. 



http : / / astronomy . swin . edu . au/ debris disks/ 
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Fig. 3. — Simulations of a system containing a single Jupiter mass planet and 500 test particles 
after 6 x 10 7 years, or 20,000 orbital periods. In each case, the planet is located at a v \ ~ 44.8 AU, 
test particles were released from parent bodies with 1.4 < a p b/a p i < 2.0, < e v \> < 0.2 and 
< i p b < 8°and /? = 0.1. The star and planet are represented by a filled star and circle respectively 
(coloured red in the online edition of the Journal). Plots (a) and (b) show the dust distributions 
generated when e v \ = 0.1 and 0.5 respectively. Plots (c) and (d) show simulated observations of 
the disks (with beam FWHM 40 AU) for e p i = 0.1 and 0.5 respectively. See the electronic edition 
of the Journal for a colour version of this figure. 
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(a) (b) 
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Fig. 4. — Simulations of a system containing a single Neptune mass planet and 500 test particles 
after 6 x 10 7 years, or 20,000 orbital periods. In each case, the planet is located at a v \ ~ 44.8 AU, 
test particles were released from parent bodies with 1.4 < a p b/a p i < 2.0, < e p 5 < 0.2 and 
< i p b < 8°and /? = 0.1. Plots (a) and (b) show the dust distributions generated when e p i = 0.1 
and 0.5 respectively. Plots (c) and (d) show simulated observations of the disks (with beam FWHM 
40 AU) for e p i = 0.1 and 0.5 respectively. See the electronic edition of the Journal for a colour 
version of this figure. 
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4.1. Dust Lifetime and Collisional Processes 

As discussed in Section 2.4, our numerical procedure of using a small number of particles 
can lead to potential problems when terminating the simulation at a specified time, as some long- 
lived particles may still be active. In a real debris disk, however, grain collisions would limit the 
maximum lifetime of particles. 

In Figure 6 we show the results of a simulation of a Jupiter mass planet with a p \ ~ 44.8 AU 
and e p i = 0.5 system containing 500 test particles with (3 = 0.1. The test particles are released from 
parent bodies with 1.4 < a^ja^ < 2.0, 0.0 < e p & < 0.2 and 0.0° < < 8.0°. The three panels show 
the dust distribution when the simulation is terminated after 11.7 million years, 40 million years 
and 153 million years, corresponding to 10%, 2% and 0% of particles remaining active. It can be 
seen that the basic disk structure is created rapidly, and the final few particles have a limited impact 
despite having lifetimes many times the median. The longest-lived particles tend to 'sharpen' the 
distribution and bring out finer detail, which is hidden in our simulated observations (telescope 
FWHM 40 AU). Thus, the imposition of a cut-off lifetime (whether as a numerical convenience 
or as a substitute for dust destruction through collisions) is of little consequence as long as most 
(>90%) particles have already been accreted or ejected. 

5. Modelling Observed Systems 

In this section we use the synthetic disk catalogue to select planetary configurations which 
generate structures similar to the observed debris disk systems Vega, e Eridani and Fomalhaut. We 
then optimise the planetary and parent body parameters, and increase the resolution (by increasing 
the number of test particles) to produce more accurate synthetic observations, which we compare 
to previous observational and numerical studies of these three debris disk systems. 

5.1. Vega 

Vega (a Lyrae) was the first star identified with an infrared excess associated with a disk 
of dusty material (Aumann et al. 1984) and has become the prototype for so-called 'Vega-type' 
stars. Subsequent observations with SCUBA (Holland et al. 1998) and the IRAM Plateau de Bure 
Interferometer (Wilner et al. 2002) confirm that Vega is surrounded by a dusty disk, with two 
prominent clumps of emission. Wilner et al. (2002) proposed that the twin lobes of emission seen 
in these images are caused by dust trapped in MMRs with a massive, eccentric planet (M p i = 
3Mj upi e p i — 0.6) orbiting Vega with a semi-major axis of 40 AU. In their model, the initial test 
particle longitudes of perihelion were constrained to be aligned closely with the planet's longitude 
of perihelion (which is reasonable for a high eccentricity system). They also noted that a wide 
variety of planetary parameters can produce the twin-lobed structure seen in Vega. 
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Whilst our code was able to reproduce the results of Wilner et al. (2002), their synthetic obser- 
vations (at the resolution of the Holland et al. (1998) Vega observations) show nearly symmetrical 
emission, whereas the Holland et al. (1998) observations show a significant asymmetry. (It should 
be noted, however, that there are uncertainties in the Holland et al. (1998) observations. In the 
model we are about to present, we are assuming that the observed asymmetry is real.) Based on 
results from our synthetic catalogue, we have modelled Vega using an entirely different planetary 
configuration in an attempt to better match these observations. In our model, a more distant 
(a p i = 73.7 AU), less eccentric (e p i = 0.1) 3 Jupiter mass planet reproduces the observed disk 
structure, with no constraints on the initial test particle perihelia (since we are modelling a lower 
eccentricity planet). We use 5000 test particles released from parent bodies with initial orbital 
elements in the range 90 < a p b < 120, 0.0 < e p b < 0.3, and 0° < i p b < 8°. 

Dent et al. (2000) state that to fit the observed spectral energy distribution, the dust grains 
which make up the Vega disk must have diameters in the range 60-400 /mi, corresponding to (3 
= 0.02 - 0.11 for Vega's estimated mass and luminosity. We have simulated this system with (3 
values of 0.02, 0.05 and 0.1, and found negligible differences between the synthetic observations. 
(Note that Wilner et al. (2002) used (3 = 0.01 with sw = 0.) Figure 7 shows the dust distribution, 
simulated observation and resonance occupancy plots obtained using our Vega model with (3 = 0.05. 

We note that our model for Vega produces a dust distribution which is essentially stationary 
in the planet's frame of reference. Such a distribution means that as the planet orbits the star, 
observations from Earth will show positional changes in the emission peaks over an orbital timescale. 
Dust distributions from the four orbital phases recorded are shown in Figure 8. 

The model accurately reproduces the twin lobes seen in the Wilner et al. (2002) observations, 
as well as the extended emission seen in the lower resolution Holland et al. (1998) observations. 
Figure 7a shows asymmetry in the two emission features, which are not collinear with the star nor 
the same distance from the star, as seen in the Wilner et al. (2002) observations. This model also 
provides reasonable constraints on the orbital parameters of the proposed planet; for example, the 
same model with e p \ — 0.2 results in a markedly different structure, as does a significantly less 
massive {M p \ < Mj up ) planet. The primary concern with this model is the requirement that such 
a massive planet have formed or migrated to such a large distance from the central star, although 
the fact that Vega is estimated to be 2.5 times as massive as the sun may mitigate this problem. 
Future detailed observations are required to test our planetary model. 

5.2. e Eridani 

First identified as a Vega-type star in 1984 (Aumann 1985), e Eridani has been the subject of 
ongoing planetary speculation since disk images taken at 850 /im by Greaves et al. (1998) showed 
a clumpy, asymmetric ring, with a cleared region interior to ~ 30 AU. Radial velocity observations 
by Hatzes et al. (2000) detected a 1.7Mj up planetary companion with a ~ 3.4 AU and e p \ ~ 0.6. 



- 17- 



Such a close companion, however, cannot account for the observed disk structure seen by Greaves 
et al. 

Numerical simulations by Ozernoy et al. (2000) have suggested the presence of a M p i = 
0.2Mj upi dpi = 55 — 65 AU, e p \ = 0.0 planet, while Quillen & Thorndike (2002) proposed a 
M p i = 0.1Mj up , dpi = 41.6 AU, e p i = 0.3 planet in the e Eridani disk. A planet with mass equal to 
or greater than Jupiter is unlikely to be responsible for the observed disk structure, since massive 
planets tend to trap particles in the n : 1 resonances (Kuchner & Holman 2003). Our simulations 
of a Jupiter mass planet, shown in Figure 9a, demonstrate this effect. For particles to occupy the 
n + 1 : n resonances near a massive planet, their parent bodies must also reside very close to the 
planet, and in this situation most particles are quickly scattered or 'leak past' the planet, as shown 
in Figure 9b. This results in a structure which lacks the observed central cavity of e Eridani (see 
Figure 9c and 9d). The presence of the observed central clearing would then require the gravita- 
tional influence of a second, more massive planet at an intermediate semi-major axis (Liou & Zook 
1999). 

Ozernoy et al. consider only particles in the 2:1 and 3:2 resonances (in equal proportions), 
while Quillen & Thorndike consider only particles with a tp = [1.1, 1.5]a p i. An examination of our 
results shows why; if particles with semi-major axes smaller than the semi-major axis of the planet 
are allowed (as in Figure 10a), emission from particles close to the star dominates and a ring 
structure is not seen (see Figure 10b). The implied ejection of particles from the inner regions of 
the disk hints at another undiscovered and more massive planet at an intermediate distance from 
the star. The requirements for such a planet are discussed below. 

We suggest that the system proposed by Ozernoy et al. is unlikely to be responsible for the 
e Eridani disk for two reasons: (1) their model is symmetrical whereas the original observations 
by Greaves et al. (1998) and subsequent 350 fim observations by D. J. Wilner (2004, private 
communication) state that only the most prominent clump (in the southeast of the image) is 
definitely present, and (2) our simulations show that dust released in a range encompassing the 
2:1 and 3:2 resonances will always occupy other resonances such as the 5:3, and that the resonance 
occupancies are rarely in equal proportion. We therefore investigate the model proposed by Quillen 
& Thorndike in more detail. 

Our models differ somewhat from those of Quillen & Thorndike in that we use 5000 test 
particles (increasing the resolution by an order of magnitude), we include solar wind drag, and 
model the system in three dimensions. The simulations were terminated after 3 x 10 8 years, when 
less than 1% of the test particles remained active. Because the dominant grain size in the e Eridani 
disk has not been well determined, we have simulated the system using two (3 values, (3 = 0.1 and 
0.01. Initial parent body semi-major axes and eccentricities are set in the range [1.1 a p i, 1.5 a p i] 
and [0,0.4] respectively, while initial inclinations are in the range [0°,8°]. Our results are shown in 
Figure 10. 

Both the (3 = 0.1 and (3 = 0.01 models reproduce the ring structure seen in the sub-mm 



- 18- 



observations, possessing a single prominent emission maxima and one or more secondary maxima. 
As expected, the contrast and detail is higher in the (3 = 0.01 case, since more particles are trapped 
in resonances for longer. Unlike Vega, this dust distribution generated by this model is not fixed in 
the planet's reference frame, but the emission peaks do show positional changes over a planetary 
orbit, as shown in Figure 11. Our results confirm the model proposed by Quillen & Thorndike 
(2002) as a viable explanation for the e Eridani system, with the caveat that an additional massive 
body must exist interior to the dust-sculpting planet's orbit. 

In order to confirm that an inner planet is responsible for clearing the inner regions of the 
disk, one would like to simply introduce a second planet into the e Eridani model. However, since 
the introduction of a second planet introduces variations into the orbital elements of the primary 
outer planet, recording particle positions with the outer planet in an fixed position over time is no 
longer possible. We can instead introduce an inner planet to our model at a range of semi-major 
axes and make use of test particle occupancy versus semi-major axis plots to determine if the inner 
planet effectively removes particles from the inner ~ 40 AU of the system. 

Plots of occupancy versus semi-major axis are shown in Figure 12 for our preferred model of 
e Eridani (Figure 12a) and models including the addition of a second, Jupiter mass planet with 
e p i = 0.3 at three different semi-major axes; a p i = 10, 15 and 18 AU (Figure 12b, 12c, and 12d). 
We see that the addition of a Jupiter mass planet with a semi-major axis between 10 — 18 AU 
is capable of substantially reducing particle occupation of semi-major axes less than that of the 
exterior planet. We found that a planet located too far from the star (a p i > 20 AU) disrupted the 
orbit of the outer planet and resulted in a dynamically unstable system. Note the change in the 
scaling of the graphs in Figure 12. We can clearly see that when an inner planet is included, the 
resonance occupancy decreases, reflecting a reduction in particle lifetimes. An understanding of the 
exact effects of a massive inner planet will have to await detailed simulations with large particle 
numbers, which can resolve the disk at all times. Such an investigation is planned for a future 
paper. 



5.3. Fomalhaut 

Unlike the previous systems, Fomalhaut has not yet been subject to detailed numerical simu- 
lations. The Fomalhaut system also differs markedly in that 850 /im SCUBA observations suggest 
the system is seen nearly edge-on and consists of a dusty torus, rather than a thin disk (Holland 
et al. 1998). The improved spatial resolution offered by the 450 fim SCUBA images of Holland 
et al. (2003) confirm the torus structure, while revealing a previously unseen arc of emission near 
or within the torus. This departure from a uniform structure is strongly suggestive of a planetary 
presence. 

The generation of a single arc of emission can be achieved by a 1:1 resonance, as suggested by 
Holland et al. (2003). Inspection of results for resonance occupancy obtained in the construction 
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of our synthetic catalogue, however, indicate that the 1:1 resonance is difficult to populate under 
normal circumstances. An example is shown in Figure 13, where a Jupiter mass planet populates 
the 1:1 resonance when parent bodies exist interior to the planet, but does not populate the 1:1 
resonance when the parent bodies are all external to the planet. We consider it unlikely that 
parent bodies with semi-major axes slightly less than the planet's semi-major axis would survive 
for the required length of time. As an alternative explanation of the single arc feature, we find 
that systems including a massive (M p i > Mj up ), moderately eccentric (0.3 < e p \ < 0.5) planet 
can induce a single emission arc over a background ring by trapping many particles in the n : 1 
resonances, where n > 1. 

After examining our synthetic catalogue, we chose a model with the following parameters: 
= 2.3M (appropriate for an A3V star such as Fomalhaut - Barrado y Navascues et al. 1997), 
M p i = 2Mj up , e p i — 0.4, and a p \ ~ 59 AU. Parent bodies were distributed with initial orbital 
parameters in the following ranges: < e p & < 0.3, 100 < a p b < 180 AU, < < 25°. The 
higher parent body inclinations naturally generates a torus structure, which is believed to exist in 
the Fomalhaut system. The simulation used 1000 test particles with (3 = 0.05, since the mass of 
dust grains in the Fomalhaut system with diameters > 100/xm (corresponding to values of f3 > 0.1) 
is believed to be negligible (Dent et al. 2000). The results of our simulations, which ran until no 
particles remained after 212 million years, are shown in Figure 14. The simulated observations bear 
a close resemblance to the 'raw' 450/xm image of Fomalhaut presented in Holland et al. (2003). 

This planetary configuration does generates a dust distribution which does not rotate with the 
planet, but appears fixed from a viewpoint external to the system over an orbital period. Thus, 
observations of Fomalhaut over time would show no change in emission if this model is correct. 
The effects of planetary phase are shown in Figure 15. 

Whilst the planetary configuration we have presented here displays a close similarity to the 
observations of Fomalhaut to date, as noted above results from our synthetic catalogue suggest that 
several parameters from the model may be varied slightly whilst still producing similar results. 
The semi-major axis of the planet is well constrained by the location of the dust ring, but the 
the production of a single arc of emission is seen with several combinations of massive, moderate 
eccentricity planets. Thus, this model is only representative of a class of planetary configurations 
which could be responsible for the structure seen in the Fomalhaut disk. 

6. Conclusions 

We have modified an N-body symplectic integrator to include the effects of Poynting- Robertson 
and solar wind drag in order to model collionless debris disks, and used the modified integrator to 
produce a synthetic disk catalogue containing approximately 300 model disks. The catalogue has 
two applications: the comparison of theoretical predictions of disk structure to numerical results, 
and to give an idea of the parameter space which might satisfy a particular observed system, 
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possibly identifying planets which are presently undetectable through Doppler shift techniques. 
Using this synthetic catalogue as a guide, we have produced 'best fit' models for three observed 
debris disk systems whose disk structures can be explained by the presence of a planetary companion 
or companions. 

We have modelled Vega using a planet with M v \ = 3Mj upi a p \ = 73.7 AU and e p \ = 0.1. The 
dust distribution generated by this model rotates with the planet, meaning future observations of 
Vega should show changes in spatial dust emission as the planetary phase changes. However, since 
a range of planetary parameters can produce the dust distribution inferred for the Vega system, 
our solution is representative of that produced by a class of massive (M p i > Mj up ), low eccentricity 
(e p i < 0.2) planetary companions. 

e Eridani was modelled using a planet with M p \ = 0.1Mj upi a p \ = 41.6 AU and e p \ = 0.3, a 
model first proposed by Quillen & Thorndike (2002). Although the dust distribution generated by 
this model does not rotate with the planet, it does change over the course of an orbital period, 
offering hope for confirmation by future observations. However, this model is dependent on the 
existence of a second, massive inner planet to clear particles interior to the modelled planet, which 
cannot be simulated using the techniques outlined in this paper. 

We modelled Fomalhaut using a planet with M p \ = 2Mj upi a p \ = 59 AU and e p \ = 0.4. The dust 
distribution generated by this model does not change significantly over the course of a planetary 
orbit, and thus the effect of planetary phase would be difficult to detect in future observations. Like 
the model presented for Vega, our model for Fomalhaut is the best match to current observations 
from a class of massive, moderate eccentricity planets which generate similar structure. 

Numerical modelling of planets in dusty debris disks has shown that planets which are unde- 
tectable through present techniques are likely to be responsible for the observed asymmetries of 
known debris disk systems. It should be noted, however, that it is difficult to positively identify a 
unique planetary companion without temporal information showing the changes in disk structure 
over the course of a planetary orbit. 

The authors wish to thanks Hal Levison for useful discussions about RMVS3 and Dave Wilner 
for providing an unpublished 350 /im e Eridani image. We also thank the anonymous referee for 
providing useful feedback and suggestions. AD was supported by a Summer Vacation Scholarship 
from the Swinburne Centre for Astrophysics and Supercomputing. All simulations were run on the 
Swinburne supercomputer 2 . 

Our synthetic debris disk catalogue is available online at: 
http : //astronomy . swin . edu . au/debrisdisks/ 



2 http : //supercomputing . swin . edu . au/ 
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Fig. 5. — Resonance occupation for a system containing a single planet with (a) M p \ = Mj upi (b) 
M p i = 0.05Mj up . In each case, the planet has e p \ — 0.1 and a p i ~ 44.8 AU. Note the lower mass 
planet tends to trap particles in narrower resonances, closer to the planet's semi-major axis. 
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Fig. 6. — Integrated dust distribution from a system containing a Jupiter mass planet with eccen- 
tricity 0.5 and semi-major axis 44.8 AU, and 500 test particles with (3 = 0.1. Parent bodies had 
initial orbital elements 1.4 < a p i,/a p i < 2.0, 0.0 < e p 5 < 0.2 and 0.0° < i p 5 < 8.0°. Simulated obser- 
vations utilised a 2D Gaussian with FWHM of 40 AU. Plots (a) and (b) show the dust distribution 
and a simulated observation after 11.7 million years (10% of particles remaining). Plots (c) and 
(d) show distribution and simulated observation after 40 million years (2% of particles remaining). 
Plots (e) and (f) show distribution and simulated observation after 153 million years (no particles 
remaining). See the electronic edition of the Journal for a colour version of this figure. 
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Fig. 7. — Results of Vega simulation using a p i = 73.7 AU, e p i = 0.1 and 5000 test particles 
with (3 = 0.05 released from parent bodies with 90 < a p b < 120 AU, 0.0 < e p 5 < 0.3, and 
0° < i v \) < 8°. The simulation was terminated after 10 8 years with 0.6% of particles remaining, 
(a) Dust distribution - note the two peaks as seen in the Wilner et al. (2002) interferometric 
observations, (b) Simulated observation of the system using a resolution of 108 AU (equal to the 
resolution of the 850 /im SCUBA observations by Holland et al. 1998) and a 5 mJy solar flux, 
with contours at 2 mJy separation. The two peaks merge together at low resolution, although one 
dominates slightly, (c) Resonance occupancy in the Vega model. See the electronic edition of the 
Journal for a colour version of this figure. 
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Fig. 8. — Dust distribution generated by the Vega simulations (with (3 = 0.05) at the four different 
recorded planetary phases. The distribution rotates with the planet, meaning observations from 
Earth will show positional changes in the emission peaks over the planetary period. 
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Fig. 9. — Structures resulting from a Jupiter mass planet with a p \ ~ 44.8 AU and e p \ = 0.3. 
Test particles have (3 = 0.1, and are released from parent bodies with initial eccentricity range 
e pb — [0.0,0.2]. (a) Resonance occupancy when the initial range of a p b is [2.0, 3.5]a p i. The higher 
order n:l resonances dominate. Note the high number of total particle positions recorded, indicating 
longer average particle lifetimes, (b) Resonance occupancy when the initial range of a p b is [1.0, 
lA]a p i. The planet rapidly removes test particles, and only a few resonances close to the planet's 
semi-major axis are populated, (c) The dust distribution created when the initial range for is 
[2.0, 3.5]a p i. A distinct ring structure is formed, containing two (unequal) density enhancements, 
(d) The dust distribution created when the initial range for is [1.0, 1.4]a p j. A ring structure is 
not generated 
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Fig. 10. — Simulations of e Eridani including a M p \ = 0.1Mj upi e p \ — 0.3 planet with a semi-major 
axis of 41.6 AU, from a viewing angle of 30°. Parent bodies had 1.1 < a p b/a p i < 1.5, < e p & < 0.3 
and 0.0 < i p b < 8.0. The telescope beam FWHM was set to 45 AU - the approximate beam 
diameter of the Greaves et al. (1998) images, (a) Dust distribution, (3 — 0.1, all particles allowed, 
(b) Simulated observation, (3 = 0.1, all particles allowed. Emission close to the star dominates and 
the ring structure is faint, (c) Dust distribution, (3 = 0.1, particles with semi-major axes less than 
the planet removed, (d) Simulated observation, /3 = 0.1, particles with semi-major axes less than 
the planet removed, (e) Dust distribution, (3 — 0.01, particles with semi-major axes less than the 
planet removed, (f) Simulated observation, (3 — 0.01, particles with semi-major axes less than the 
planet removed. See the electronic edition of the Journal for a colour version of this figure. 
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Fig. 11. — Dust distribution and synthetic observations from by the e Eridani simulations (with 
(3 — 0.1 and particles with semi-major axes interior to the planet removed) at three different 
planetary phases. Dust distributions are shown in (a) to (c), with the corresponding synthetic 
observations shown in (d) to (f). The distribution does not rotate with the planet, but the emission 
peaks show positional changes over the course of an orbit. See the electronic edition of the Journal 
for a colour version of this figure. 
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Fig. 12. — Occupancy versus semi-major axis plots showing the effect of adding a Jupiter mass, 
e p i = 0.3 inner planet to the e Eridani model. Plot (a) shows the plot for the unmodified e Eridani 
model - note test particles occupancy of semi-major axes interior to 40 AU. Plots (b), (c) and 
(d) show the results for systems with an inner planet added with a semi-major axis of 10, 15 
and 18 AU respectively. In each case, test particle occupancy of semi-major axes interior to 40 
AU is substantially reduced (note the different scales on the four graphs). Simulations contained 
five hundred (3 = 0.1 test particles were released from parent bodies with 46 < a p b < 62 AU, 
0.0 < e v \) < 0.4 and 0.0° < i v \) < 8.0° and were terminated after 10 8 years. 
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Fig. 13. — The resonance occupancy for a 500 particle simulation of a system containing a Jupiter 
mass planet with a v \ = 40 AU and e v \ = 0.1. Test particles had (3 = 0.1. (a) With parent bodies in 
the range 35 < a p & < 45 AU, 0.0 < < 0.3 and 0.0° < i v \) < 8.0°, the 1:1 resonance is prominent, 
(b) When parent bodies are in the range 45 < a p b < 60 AU, 0.0 < e p 5 < 0.3 and 0.0° < i p b < 8.0°, 
there is no population of the 1:1 resonance. 
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Fig. 14. — Simulations of the Fomalhaut system using a 2 Mj up planet with e p i = 0.4 and a v \ ~ 59 
AU. (a) The resonance occupation - note that the n:l resonances dominate, (b) Face-on distribution 
of particles, (c) Simulated observation of system at 450/xm from an inclination of 70°. The telescope 
FWHM of is 50 AU, equivalent to the resolution of SCUBA at 450/im. Contours show equal intensity 
increments. See the electronic edition of the Journal for a colour version of this figure. 
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Fig. 15. — Dust distribution generated by the Fomalhaut simulations at three different planetary 
phases. The distribution does not rotate with the planet, and remains relatively constant in a fixed 
reference frame. 



